%% 1. The Hansen model as the basline comparsion (different from paper model basline)
par = par_ini;

disp('-------------> Model 1, Hansen (1985) <----------------------')
disp('Note: adjusted by aggregate-level of productivity, taxes and G<----------------------')

B_Hansen       = (1 / par.BETA - (1 - par.DELTA)) / (1 - par.tau_k) / par.q_average_prod;
w_Hansen       = par.ETA / (B_Hansen / par.A / (1 - par.ETA))^((1 - par.ETA) / par.ETA);
c_star_Hansen  = (par.XI / (1 - par.tau_h) / w_Hansen)^(1 / (-par.SIGMA));
K_Hansen       = (c_star_Hansen + par.G) / (B_Hansen * par.q_average_prod / (1 - par.ETA) - par.DELTA);
K_bar_K_Hansen = par.q_average_prod;
Y_Hansen       = B_Hansen * par.q_average_prod * K_Hansen / (1 - par.ETA);
H_Hansen       = (Y_Hansen / (par.A * K_bar_K_Hansen * K_Hansen)^(1 - par.ETA))^(1 / par.ETA);
Util           = log(c_star_Hansen) - par.XI * H_Hansen;
vec            = [B_Hansen, K_Hansen, K_bar_K_Hansen, c_star_Hansen, H_Hansen, Y_Hansen, 0, 0];
C              = {'B','K','K prod','c','H','Y','AQC','SPPE'};
D              = C;
D(2,:)         = num2cell(vec);
disp('The steady state result is')
fprintf('%10s %d\n',D{:})
disp(' ')

%% 2. The Hansen model + idiosyncratic risks

%%%% Notice that we have mean preserving spread on the productivity so the
%%%% aggregate outcome is the same as in the Hansen model above (check the
%%%% result shown in the Command Window)

disp('-------------> Model 2, Hansen + idiosyncratic shocks but no reallocation<----------------------')
disp('Mean preserving spread; aggregate allocation is the same as Hansen Model')

par           = par_ini;
par.ALPHA     = 0;
fun           = define_function_fun(par);

perfect_reallocation = 0;
cali_step     = 0;
x             = fsolve(@(x)eqm_ss_iid_cases_2_to_5(x, par, fun, perfect_reallocation, cali_step), [0.15, 0]);
[fval, eqm]   = eqm_ss_iid_cases_2_to_5(x, par, fun, perfect_reallocation, cali_step);
B             = x(1);
ALPHA_0       = x(2);
vec           = [B, eqm.K, eqm.K_bar_K, eqm.c_star, eqm.H, eqm.Y, 0, 0];
C             = {'B','K','K prod','c','H','Y','AQC','SPPE'};
D             = C;
D(2,:)        = num2cell(vec);
disp('The steady state result is')
fprintf('%10s %d\n',D{:})
disp(' ')
prod_std_Hansen = sqrt( (exp(par.log_eps_sig^2) - 1) * exp(2 * par.log_eps_mu + par.log_eps_sig^2) ); %log normal distribution's STD.DEV. 